Multiplicación de Matrices con un kernel en bloques: una Introducción a programación con teselas (Tiled Programming)

En esta sección nos dedicaremos a introducir un método de programación mucho más eficiente. En inglés es llamada tiled programming y no hay una traducción exacta en español. Por eso la llamaremos programación con teselas. Estas últimas vienen a ser la traducción de Tile, que es un tipo de bloque. Y puesto que en CUDA ya tenemos los bloques de threads o simplemente bloques, el nombre de programación en bloques viene a ser confusa y ambigua, puesto que las teselas, a diferencia de los bloques, contienen información (datos) y no threads.

Por lo que de ahora en adelante habrá de recodar que cuando se habla de programación con teselas, estas últimas se refieren a bloques de información y no de threads.

En el ejemplo concreto del anterior kernel de multiplicación de matrices, cada dato de las dos matrices se encuentra en la memoria global y es utilizado por el kernel en distintas ocasiones. La idea de la programación con teselas es reducir y ordenar el número de accesos a los datos con ayuda de los distintos tipos de memoria que se encuentran en el GPU. Esto se observa mejor en las figuras siguientes.

La ilustración superior muestra el proceso cuando no programamos agrupando información. En las imágenes inferiores pasamos a un proceso por `teselas`.

En estas dos imágenes observamos como se hace el cálculo usando teselas. El transporte de información se observa más limpio y ordenado

Supongamos entonces que los threads 1 y 2 utilizan 6 datos guardados en la memoria global. Para optimizar los cálculos, estos los haremos en 2 fases. En la primera guardamos sólo 3 datos en un bloque de la memoria compartida. Ahora los threads 1 y 2 tendrán un acceso más rápido a los datos que si estos estuvieran en la memoria global. Sin embargo hay que recordar que la memoria compartida tiene un tiempo de vida igual al thread block, por lo que cuando los threads 1 y 2 terminen con los 3 datos, estos serán eliminados de la memoria compartida. Ahora entramos en la fase 2, en la que guardamos los siguientes 3 datos en la memoria compartida. repitiéndose lo mismo que en la fase 1.

Siguiendo la analogía...

Pensemos ahora que cada uno de los niños ya está en su casa después de un viaje en el camión escolar y están trabajando arduamente en sus tareas. Sin embargo ahora son sus padres quienes se preocuparán búscando qué preparar para cenar, por lo que habrán de ir al super mercado para comprar los ingredientes con los que se preparará la comida.

En muchas ciudades del mundo, la distribución de algunos productos se encuentra centralizado en algunos barrios de la ciudad, por lo que encontrar algunos productos requiere de viajar grandes distancias hasta el lugar en el que se encuentra el producto deseado. Sin embargo desde hace décadas algunas otras ciudades han implementado políticas para descentralizar estos mercados, por lo que si una familia vive en algún barrio suburbano, esta no tenga que viajar hasta el centro de la ciudad para conseguir algunos productos. Esto le permite a las personas optimizar sus viajes.

La idea de la programación en bloques es la misma. Los threads necesitan operandos (comida) para trabajar. Al programar en bloque uno lleva la comida al super mercado local (memoria compartida) en el que se encuentra el thread de tal manera que este no tenga que viajar al centro de la ciudad (memoria global).

Sin embargo podemos encontrar dos distintos tipos de problemas que pueden presentarse en estos casos. Uno es cuando en el barrio local, cada una de las familias que viven ahí no buscan los mismos productos, es decir que en un thread block, cada uno de estos no busca los mismos operandos. Otro es cuando el supermercado local no tiene la suficiente infraestructura para darle cábida a toda la comida necesaria, i.e. no hay el suficiente espacio en la memoria compartida para traer todos los operandos necesarios.

En el primer caso, hay que encontrar una serie de datos la cual sea utilizada por múltiples threads, asegurándonos así que cada dato sea utilizado, para luego pasar a la siguiente serie de datos una vez que los threads hayan hecho sus operaciones.

El segundo caso tiene que ver con la capacidad de la memoria compartida. Podemos darle a esta dos configuraciones: 16 KB de memoria compartida o 48KB en el cache L1. Esto nos permite dar de nuevo los detalles a la hora de elegir el tamaño del bloque, o viceversa. Tomemos aquella con 16KB de memoria compartida.

  • Supongamos que elegimos de nuevo un tamaño de thread block de 16x16. Puesto que cada thread trabaja con dos datos (uno de la matriz $A$ y otro de la matriz $B$) y cada uno de estos pesa 4B, entonces obtenemos que cada uno de los bloques utiliza 256x2x4B = 2KB, por lo que podremos utilizar los 8 bloques permitidos para llenar los espacios permitidos en el SM.
  • Si ahora utilizamos un bloque de 32x32, entonces cada bloque usará 1024x2x4B = 8KB, resultando en que únicamente podremos utilizar 2 bloques. Pero como ya habíamos visto en el notebook anterior, al utilizar una dimensión de 32x32 entonces estaremos limitados a usar un único bloque.

De esta manera la dimensión del bloque y el número de operaciones hechas por los threads son a considerarse a la hora de utilizar la programación en bloques.

Después de todas estas restricciones y detalles en los que hay que pensar, uno llega a preguntarse: ¿realmente vale la pena este tipo de programación?

Por ahora no tenemos las herramientas suficientes para responder la respuesta. Sin embargo, esperamos que el lector tenga confianza y no piense en que los autores escribieron este notebook sólo para hacer perder el tiempo.

Por ahora nos limitaremos a re-escribir el código de multiplicación de matrices usando el método con teselas. Una vez que tengamos los dos códigos entonces podremos comparar dos kernels que realizan la misma función con dos métodos distintos.

Multiplicación de matrices usando teselas

La idea general del código es bastante sencilla. Supongamos que buscamos multiplicar dos matrices $A$ y $B$ de 16x16. Sin el método con teselas la multiplicación se daría de manera simultanea con todos los datos de las dos matrices, tal cual se muestra en la primera figura de este notebook.

Con el método con teselas, cada matriz se puede dividir -digamos- en 16 teselas de 4x4. De esta manera tomaremos una tesela de la matriz $A$ y otra de la matriz $B$ para obtener aquella porción de la matriz resultante que puede ser calculada con las dos teselas. De esta manera, la combinación de teselas de las dos matrices van permutando hasta finalmente sumar la matriz resultante completa. Todo esto hecho simultaneamente. La imagen siguiente da muestra de esto.

Esquema en el cual se muestra como interactúa una tesela de la matriz A con dos teselas de la matriz B para calcular una parte de la matriz C.

Sin más, comencemos.

Lo primero que tenemos que hacer en nuestro código será declarar las teselas (o tiles) en la memoria compartida. Esto se hace de la siguiente manera:

#define ANCHO_TESELA 16

__shared__ float tesela_A[ANCHO_TESELA][ANCHO_TESELA] ;
__shared__ float tesela_B[ANCHO_TESELA][ANCHO_TESELA] ;

Se define ANCHO_TESELA afuera del kernel, pero aquí lo ponemos para fines demostrativos. Lo que hace ANCHO_TESELA es fijar la dimensión de nuestros bloques de información, que en este caso será 16 por discusiones anteriores (nb 5 y unos párrafos atrás) sobre la eficiencia de nuestro código.

Para llenar la memoria compartida de los datos de nuestra matriz, tendremos que usar un bucle for con el que podamos pasar por todos los bloques de threads que se tienen.

for (int t = 0 ; t < (numColA - 1)/ANCHO_TESELA + 1 ; t ++) 
{

}

donde t es el número de bloque. Es importante notar la expresión (numColA - 1)/ANCHO_TESELA + 1, que como resultado nos da el número de separaciones que la matriz puede tener con la dimensión de las teselas. Dicho de otro modo, el parámetro t utilizado en el for nos da el número de tesela en el cual estamos trabajando. El número máximo de tiles con los que trabajaremos está dado, en el caso de una multiplicación de matrices de $m\times n$ y $n\times k$ por (n-1)/ANCHO_TESELA + 1; algo así como el algoritmo que utilizamos en el capítulo anterior para saber la dimensión de nuestra grid.

Nota: Esto nos da una buena idea de las dimensiones de la malla y de los bloques que luego se tendrán que fijar en el código.

Dentro de este for habrá que llenar la memoria compartida. Esto se hace con un código de la forma

if (Fila < numFilasA && t*ANCHO_TESELA + tx < numColA) 
{
    tesela_A[ty][tx] = A[Fila*numColA + t*ANCHO_TESELA +tx] ;
} else {
    tesela_A[ty][tx] = 0.0 ;
}

Los dos primeros if/else son condiciones de frontera y nos permiten almacenar la cantidad exacta de datos. Supongamos que nuestra matriz $A$ es de 256x257, entonces necesitaremos de 17 teselas de dimensión 16 en la dirección x, haciendo que 15 threads no tengan asignados ningún valor de la matriz $A$. Para evitar que tengan información no deseada, hacemos a todos estos sobrantes 0.0. Pasa lo mismo para la matriz $B$. El índice utilizado en la matriz $A$ es un poco más complicado que todos aquellos índices vistos hasta ahora, sin embargo le dejaremos al lector su comprensión.

En el caso de tesela_B, el código es el siguiente:

if (t*ANCHO_TESELA + ty < numColA && Col < numColB)
{
    tesela_B[ty][tx] = B[(t*ANCHO_TESELA +ty)*numColB + Col] ;
} else {
    tesela_B[ty][tx] = 0.0 ;
}

Recuerda, a la hora de escribir el código, de colocar un __syncthreads() al final de estos dos párrafos para evitar cálculos equivocados.

Finalmente calculamos las entradas de la nueva matriz.

for (int i=0; i < ANCHO_TESELA; i ++) 
{
    c_resultado += tesela_A[ty][i] * tesela_B [i][tx] ;
}

__syncthreads() ;
if (Fila<numFilasA && Col<numColB) 
{
    C[Fila*numColB + Col] = c_resultado ;
}

En esta pequeña parte del código se calculan las entradas de la matriz resultante $C$. Nota como se utilizan dos variables, C y c_resultado. La primera variable es el arreglo correspondiente a la matriz $C$, mientras que c_resultado es un número flotante correspondiente a las entradas de la matriz $C$.

Observando con atención la esencia del código de multiplicación de matrices anterior no ha desaparecido, y únicamente agregamos unos pasos intermedios para optimizar nuestro código.

Ahora dejamos el lector para que pueda experimentar con estos nuevos conocimientos. Como ejercicio tendrá que escribir un código en el cual se realice la multiplicación de dos matrices con el método con teselas. Las matrices son las mismas que en el notebook pasado para poder comparar los tiempos y resultados.

Sin más, ¡a escribir se ha dicho!


In [ ]:
%%writefile Programas/Mul_matrices_teselas.cu

#define ANCHO_TESELA

// Escribir el kernel

// En primer lugar tendras que declarar las teselas
// Luego copiar los datos de la memoria global a cada una de las teselas
// Para luego poder calcular los datos de la matriz resultante

int main( int argc, char * argv[])
{
    
    int numColA 1000 ;
    int numColB 1000 ;
    int numColC 1000 ;
    
    int numFilaA 1000 ;
    int numFilaB 1000 ;
    int numFilaC 1000 ;
    
    float A[numFilaA*numColA] ;
    float B[numFilaB*numColB] ;
    float C[numFilaC*numColC] ;
    
    for (i = 0 ; i < numFilaA*numColA ; i ++)
    {
        A[i] = i ;
    }

    for (i = 0 ; i < numFilaB*numColB ; i ++)
    {
        B[i] = i + 1 ;
    }
    
    // Escribe abajo las lineas para la alocacion de memoria
    
    // Escribe abajo las lineas para la copia de memoria del CPU al GPU
    
    // Escribe las dimensiones de los bloques y mallas.
    // En el caso de la malla, la dimension esta dada en este caso por (n-1)/ANCHO_TESELA + 1
    
    // Escribe las lineas para lanzar el kernel
    
    // Copia la memoria del GPU al CPU
      
    // Libera la memoria
      
    // Escribe un codigo para poder saber si tu resultado es correcto
    
    
    return 0 ;
}

In [ ]:
!nvcc -o Programas/Mul_matrices_teselas Programas/Mul_matrices_teselas.cu

Probado localmente, se utilizaron dos matrices con dimensiones 200x100 y 100x256. Al correr nuestro código con el kernel básico, el cálculo se realizó en 0.24 ms. Cuando se compiló el código con un kernel en bloques o tiled kernel, el computo se realizó en ¡¡¡0.13ms!!!.

Ejercicios

[1] Implementar funciones para comparar los tiempos entre los dos kernels y probar la eficacia de la programación con teselas.

[2] Jugar con las dimensiones de las teselas para mostrar cómo cambia la eficacia del código con respecto a este parámetro.

Conclusiones generales de CUDA C

En esta primera parte nos hemos apoyado de ejemplos clásicos para poder mostrar las bases de la programación en CUDA C. A partir de aquí quedan todavía una enorme cantidad de detalles por ver, sin embargo consideramos que el lector, para los fines de este curso, podrá continuar sin mayor problema.

En el siguiente notebook hemos preparado una serie de ejercicios para poder seguir con el entrenamiento de CUDA C. También, para ir más lejos en algunos tópicos invitamos a consultar los libros puestos en las referencias.

CUDA C, aunque parezca tedioso, tiene un gran número de ventajas. Esperamos que hasta este momento el lector se haya ya convencido de todo aquello que se puede hacer con la computación en paralelo y que decida continuar con nosotros para ahora empezar con cosas más interesantes.

En la segunda parte usaremos la física estadística como pretexto para aprender aún más de CUDA C y de pyCUDA, e inversamente nos dedicaremos a enseñar un poco de física estadística usando las ventajas de CUDA C y pyCUDA.


In [ ]: